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J5 ■ ABSTRACT 

Context. Ideal magnetohydrodynamical (MHD) simulations have revealed catastrophic magnetic braking in the protostellar phase, 

■ which prevents the formation of a centrifugal disk around a nascent protostar. 
Aims. We determine if non-ideal MHD, including the effects of ambipolar diffusion and Ohmic dissipation determined from a detailed 

, _ chemical network model, will allow for disk formation at the earliest stages of star formation. 

, Methods. We employ the axisymmetric thin-disk approximation in order to resolve a dynamic range of 9 orders of magnitude in 

^0 ■ length and 16 orders of magnitude in density, while also calculating partial ionization using up to 19 species in a detailed chemical 

equilibrium model. Magnetic braking is applied to the rotation using a steady-state approximation, and a barotropic relation is used 
to capture the thermal evolution. 

Results. We resolve the formation of the first and second cores, with expansion waves at the periphery of each, a magnetic diffusion 
^ ■ shock, and prestellar infall profiles at larger radii. Power-law profiles in each region can be understood analytically. After the formation 

■ of the second core, the centrifugal support rises rapidly and a low-mass disk of radius a 10i?o is formed at the earliest stage of star 
I formation, when the second core has mass ~ 10"'' Mq. The mass-to-flux ratio is ~ 10'* times the critical value in the central region. 

, Conclusions. A small centrifugal disk can form in the earliest stage of star formation, due to a shut-off of magnetic braking caused 

by magnetic field dissipation in the first core region. There is enough angular momentum loss to allow the second collapse to occur 
directly, and a low-mass stellar core to form with a surrounding disk. The disk mass and size will depend upon how the angular 
(N- momentum transport mechanisms within the disk can keep up with mass infall onto the disk. Accounting only for direct infall, we 

^ , estimate that the disk will remain < 10 AU, undetectable even by ALMA, for a 4 x 10* yr, representing the early Class phase. 

, Key words, magnetohydrodynamics (MHD) - protoplanetary disks - stars: formation - stars: magnetic field - accretion, accretion 
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1. Introduction infall ing gas. This is the classical angular momentum problem 

("e.g.. lSpitzerlll978h . The required reduction of angular momen- 
Near the beginning of the star formation process, prestellar j, ^ften credited to magnetic braking, which acts during the 

cores are observed to be rotating (e.g., — ° * 



[1981 [Goodman et al.l[T99l ICaselU et al 



Goldsmith & Arquilla| contraction and collapse by linking the core with its envelope 

| 2002|) . At the end of transfprrind andnlar mnmpntnm (<ip.p. iRasii Mniisrhnviasl 



. ^ , and tr ansferring angular momentum (see iBasu & MouschoviasI 

. the process, observations show nearly all Class II objects to references therein). Recent numerical simulations of 

^ . be surrounded by rotating disks of size > 100 AU, likely in protostellar collapse under the assumption of ideal magneto- 

; centrifugal balance undergomg Keplenan rotation (see, e.g., hydrodynamics (MHD) have suggested the converse problem, 

■ - - [ Andrews & Williams || 2005|: [ Williams & Cieza| [ 201 U and refer- „^^gjy ^^at cores may experience catastrophic magnetic brak- 

ences therein). Observational data concernmg what happens in ^^is is the result that an extremely pinched magnetic field 

between those two snapshots is relatively scarce. Currently, there configuration has enough strength and exerts a long enough lever 

IS no evidence for the presence of centrifugal disks larger than j^e inner regions of collapse that magnetic braking can 

^JO AU around Class or Class I objects (e.g., [ Maury et al. | suppress the formation of a centrifugal disk entirely. 
120101) . However, there are outflows observed around these young 

objects, commonly linked to disk accretion. Disks must therefore Catastrophic m agnetic braking was first demonstrated by 

form with a small size and grow over time. The dearth of obser- lAllen et al.) (l2003l) . who used two-dimensional (2D) axisym- 

vational data at smaller scales will be remedied by ALMA. In metric calculations . Subseq uent ideal MHD simulations by 

the meantime, the evolution of angular momentum and the for- iMellon & 11 (|2008|) in 2D and lHennebelle & Fromang|(l2008h in 

mation of a centrifugal disk have to be studied theoretically and three dimensions (3D) showed that catastrophic magnetic brak- 

numerically. ing occurs for initially aligned (magnetic field parallel to rota- 

In order to form stars of the observed sizes and rotation rates, tion axis) rotators in which the magnetic field is strong enough 

most of the angular momentum has to be removed from the that yu < 10, where yu is the initial core mass-to-flux ratio nor- 
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malized to its critical value for gravitational collapse. Disk for- 
mation in ideal MHD then requires at least ju > 10, which is 
significantly greater than the typica l value ~ 2 for dense core s 
determined fr om observations (e.g .. FlYoland & Crutcheill2008h . 
Furthermore, iMellon & Lil m6% found that the inclusion of 
non-ideal effects of Ohmic dissipation and ambipolar diffusion 
in a 2D model was not sufficient to enable disk formation on 
scales of > 10 AU from the center that were resolvable in their 
simulation. In a newer 2D mode l that also includes the Hall 
term of non- ideal MHD, iLi et al.l (1201 ll) again find essentially 
the same result. All of the above simulations have an inner res- 
olution limit of a; 10 AU due to the use of a sink cell or an 
artificially stiff equation of state at high densities. Since these 
simulations also terminated soon after the central stellar core 
formed, the main common result of these simulations is that 
large ~ 100 AU disks (as observed around Class II objects) can- 
no t form during the early Class phase. A 3D ideal MHD model 
bv lHennebelle & Ciardil(l2009l) that starts with a non-aligned ro- 
tator, also finds catastrophic magnetic braking on these scales 
but within a smaller range of fi < 3. The effect of numerical 
reconnection may account for some of the differences amongst 
these models. All of the above works leave open the possibility 
that a small disk may actually form within the inner ~ 10 AU 
during the earliest phases of star formation. 

An interesting question is: can a very small disk form in the 
innermost regions (near the protostar) in the case of ideal MHD? 
iGalli et all (l2006h developed a simplified analytic model with 
a split-monopole magnetic field to show that catastrophic mag- 
netic braking extends to the s t ellar s urfac e in the flux-frozen 
limit. iKrasnopolskv & Konigll (|2002|) and [Braiding & Wardld 
(|20I ih employed self-similar non-ideal MHD models in the 
thin-disk approximation, using t he prescription for steady-stat e 
magnetic braking developed by iBasu & MouschoviasI (Il994h . 
Depending on different levels of effective diffusivity, they 
found either catastrophically -braked or disk-formation solutions. 
However, the diffusivity was parametrized in a way that pre- 
served self-similarity, and its magnitude and functional form did 
not make contact with microphysical models. Consequently, a 
prediction for whether or not disks may actually form was not 
possible. 

Our work fills the gap of numerically studying disk forma- 
tion in the early Class phase do wn to scales smaller than a 
stellar radius. In a previous paper (lDapp&Basull2OI0l here- 
after DBIO) we followed the evolution of a collapsing molecu- 
lar cloud core in axisymmetric thin-disk geometry all the way to 
protostellar densities (« > 10^*^ cm"^). We found that Ohmic dis- 
sipation alone (in a simplified form) reduces the magnetic field 
strength sufficiently to render magnetic braking within the first 
core ineffective. As a result, a centrifugal disk can indeed form 
around the second core (the protostar). 

The only other related MHD models that res olve the stellar 
core are the 3D non-ideal MHD simulat i ons by iMachida et alJ 
(I2006ll2007h and lMachida & Matsumotol (1201 lb . who also used 
a simplified treatment of Ohmic dissipation. In most of their pub- 
lished models, the evolution was halted at up to ten days after 
the second core was formed, due to time-step restrictions, and 
the fo rmation of a centrifugal disk w as not studied. Very recent 
work dMachida & Matsumotol I20T1I) follows disk formation in 
one MHD model, and its relation to our work is discussed in 
Section [8] 

In this paper, we increase realism by employing a de- 
tailed chemical network including the effect of grains to 
calculate partial ionization and the resulting Ohmic and 
ambipolar diffusivities. We use the method presented in 



iKunz & MouschoviasI (l2009l) . Both Ohmic dissipation and am- 
bipolar diffusion gain importance with the formation of the 
first core (e.g.. |Li & M cKee 1996; Contopoulos et al. 1 1 9981 : 
iDesch & MouschoviasI I2OOII) . We follow the evolution of the 
core to stellar sizes, using a multi-fluid approach that includes 
important inelastic collisions between gas-phase species and 
grains. We also study the effect of different grain sizes. We model 
four gas-phase species and up to five different grain sizes with 
three charge states (singly positively- or negatively charged, and 
neutral) for each grain size. We demonstrate that Ohmic dissipa- 
tion and ambipolar diffusion can reduce the efficiency of mag- 
netic braking in the inner regions enough to allow for the forma- 
tion of a centrifugally-supported disk. There is dramatic mag- 
netic flux loss, and a small disk can form in a very early phase 
of evolution. We propose a disk formation scenario in which 
the disk remains < 10 AU during the Class phase, and sub- 
sequently expands significantly due to internal angular momen- 
tum transport mechanisms such as m agnetorotational instabil- 
ity (MRJ) or gravitation al torques (e.g.. lBalbus & Hawlevlll99ll : 
IVorobvov & iBasull2007h . 

This paper is structured the following way. In Section |2] 
we formulate the problem and describe the method of solution. 
Section [3] discusses the implementation of magnetic braking, 
while Section|4]outlines the derivation of the induction equation. 
Section|5]details the chemical model used to calculate the abun- 
dances of each species and Section|6]contains the description of 
our initial conditions. Readers mainly interested in our results 
may proceed directly to Section |2l The discussion is found in 
Section[8] and we summarize our findings in Section|9] 



2. Method 



We solve the equations of non-ideal MHD for a rotat- 
ing, axisy mmetric thin disk made up of partially-ionized gas 
(iBasu & M ouschovias 1994). In the thin-disk approximation, 
all quantities are assumed to be uniform over one scale height 
of the disk. It is applicable as long as the half-thickness of 
the disk remains small compared to the radius of the disk. 
iFiedler & MouschoviasI ( Il993h showed that an initially cylin- 
drical cloud quickly flattens parallel to the magnetic field be- 
fore any significant collapse in the radial direction occurs, and 
a oblate cloud is forme d that can be described by t he thin-disk 
approximation (see also lBasu & Mouschovias|[l994l) . While this 
approximation breaks down at the interfaces of the first core and 
in the second core, it holds well in most of the cloud. It allows 
us to follow the evolution to very high densities and very small 
scales with realistic physics (non-ideal MHD and sophisticated 
chemical model), and does not require a lot of computation time. 

We solve the following system of equations, derived by in- 
tegrating the MHD equations over the z-direction and assuming 
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axisymmetry: 
dt 

dt 

dL 
dt 

dt 



r or 
1 d 

: - - — (rSn V,V,) + /p + /g + /m + f„ 

r dr 

Id, , 1 

' —irLvr) + —rB^^eaB^,z, 

r dr 2n 

15/ X 
1 d I dB, 



r dr \ dr 



(la) 

(lb) 

(Ic) 

(Id) 

(le) 



Here, is the mass column density of the neutrals, iv is 
the radial velocity, and L = S„r2r- is the angular momentum 
per unit area. The mass volume density is denoted by g^. The 
forces /p, /g, /m, and /r, as well as the magnetic field compo- 
nents (B;,eq, fif.z, and B^p j) are discussed below (see Eqs.© 
as well as ( iTal i and (fTbb). We do not include a separate con- 
tinuity equation for the grains since no significant change of 
the dust-to -gas ratio occurs beyond a number density of » ^ 
lOfcm"^ (ICiolek & Mouschoviaslll996l: iKunz & MouschoviasI 
I20T0I) . Due to the very low ionization fraction in the modeled 
core, we neglect the inertia of all species but the neutral molec- 
ular gas. 

We include the ef f ect of both ambipo lar diffusion (e.g., 
iMestel & Spitzej|1956l: iMouschoviasI [19791) and Ohmic dissi- 



pation (e.g.. iNakano et al 



in an effective resistivity ?7eff, 
which we derive from microphysical considerations in Section 
m Note that we do not assume the resistivity to be spatially con- 
stant by p ulling it outside all sp atial derivatives, which is dif- 
ferent from lMachida et all (l2007h . We ignore the contrib ution of 
the th ird non-ideal MHD eff'ect — the Hall effect (e.g.. lWardl3 
l2007h — in this Paper, and leave that for future work. 

For computational ease, we use the barotropic relation 
Eq.lfTeb instead of a detailed energy equation. 

Foll owing DB 10, we parametrize t he temperature-density re- 
lation of iMasunaga & Inutsukal (l2000l) in the following way: 



i 10 



n\ < n < n2, 

^^{T^riitr «2<«<«3, (2) 
iio(r(ir(r(r 

In this expression, the temperature T is given in units of 
K. The number densitiefl « at which the powerlaws in the 
temperature-density relation change are «i = 3 x 10"^ cm"-', 
n2 = 10'^^ cm-\ «3 = 3 X 10'^ cm-\ and «4 = 10^' cm^l 

We transform this temperature-density relation into a 
pressure-density relation using the ideal gas law 



P = nkjiT 



(3) 



Here, P is the thermal pressure and is Boltzmann's constant. 
We calculate the thermal midplane pressure of the neutrals in our 



' Note that during hydrogen dissociation the mass per particle is re- 
duced from 2.3 ;«h (which is calculated assuming 10% He in number) 
to 1.3 mn, and then finally to 0.6 ran after the gas is fully ionized; is 
the mass of a hydrogen atom. 



thin disk self-consistently, including the effects of the weight of 
the gas column, external pressure, magnetic pressure, and the 
effect of a central star (once present, with mass M^): 



2 n 

+ ^ext + -3— + 



Stt 



(4) 



Here, W^^ is the extra vertical squeezing due to the star's gravity 
integrated over the disk's finite thickness. The external pressure 
is fixed at Pext =0.1 ttCLq/I, where G is the gravitational con- 
stant and So is the initial central neutral column density. We set 
- 0, and interpolate gn from the pressure-density relation 
found earlier. If a central object is present (after the formation 
of the second core, and after introduction of the sink cell), we 
then insert the so-determined approximation to g„ and the disk's 
half-thickness Z = Sn/2pn into 



W^^2 GMi,g 



Jo 



(r2+z2) 



2^3/2 ' 



(5) 



where z is the vertical coordinate, and iterate this procedure un- 
til converges. The result also allows us to find the midplane 
temperature of our disk using Eq.Q. It is very similar to that 
found in an axisym metric (r, z) simulation with explicit radia- 
tive transfer added (iKunz & Mouschoviasll2010l) . Note that the 
extra squeezing from the protostar also flattens the disk further, 
helping to justify the thin-disk approximation. 

Once the the second core forms, at a central number den- 
sity of ^ lO^" cm ^, we introduce a sink cell of size ^ i Rq 
(which is only slightly larger than the second core). We remap 
the computational domain to a new grid with the sink as the in- 
nermost cell. The first neighboring cell outside the sink has an 
extent of only » 10"^ AU ~ 0.2 Rq, and each subsequent cell in- 
creases in size by a constant factor. We distribute matter within 
the sink cell so that much of it goes into a central point mass, 
but with enough remaining in the cell so that its column density 
remains equal to that in the first neighboring cell (because the 
gravitational solver requires non-zero mass in the first cell). We 
also impose no-outflow boundary conditions on the sink cell. We 
modify the gravitational potential by adding the contribution of 
the central point mass (with mass M^). 

The forces per unit area appearing in the momentum equa- 
tion Eq.lfTbli are given by: 



dB^ 



dr 

/g = ^ngr 
/r 



(6a) 
(6b) 



(6c) 
(6d) 



These expressions are derived in lCiolek & Mouschoviasl(ll993h . 
Equation ( l6al ) is the thermal force in the thin-disk formulation 
due to pressure gradients, while Eq.(l6b]i describes the radial 
gravitational force. Equation ( l6c] i represents the magnetic force 
due to azimuthal and poloidal components of the magnetic field 
at the top and bottom surfaces of the cloud, as well as from varia- 
tionsin scale height. This expression is derived in Ciolek & Basul 
( l2006l) . Finally, Eq.(|6dli is the centrifugal force in an axisym- 
metric thin-disk geometry, and couples the force equation with 
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the equation for the evolution of the cloud's angular momentum, 
Eg. ([Tell. Note that the effects of external pressure, magnetic pres- 
sure, and squeezing due to the protostar enter the pressure force 
through the half-thickness Z, and not explicitly in Eq.(l6al). 

We calculate the gravitational acceleration gr with the 
integral method for infinites imall y-thin disks e mploy ed in 
ICiolek & MouschoviasI (Il993h and Morton et al.l (Il994l) . This 
produces a diverging field at the disk edge, which we avoid by 
adding in the gravitational effect of a virtual column density pro- 
file oc r ' extending to infinity. We also correct for the finite 
thickness of the flattened core so as not t o overestimate the field 
strength (see lBasu & Mouschovias|[l994l) . 

The magnetic field points solely in the vertical direction in- 
side the disk, and possesses radial and azimuthal components 
(Br,z and B^ z) at the disk surfaces (indicated by subscript Z) 
and beyond. B^.z is determined from a potential field assuming 
force-free and current-free conditions in the external medium, 
using the same integral kernel A4(r, r') as for the gravitational 
field. We calculate Zj^ z and implement magnetic braking as in 
iBasu & Mouschoviaa ( ll994l) (for details see Section[3]). 



The remaining constituting equations are: 

dr'r' [B,,eq(r')-Bref]M(n/). 

- V^^eext o (r) 



■[Q(r)-Q,ef], 



<D(r)= r drV'B,,eq 
Jo 

/-*oo 

grir)^27TG drV%(r')A^(r,r'), 
Jo 



Mir,r') 
Z{r) 



2dl^^ 
n ar r> \ r> 



2gn (r) 



(7a) 

(7b) 
(7c) 
(7d) 
(7e) 
(7f) 



Here, B,ef is the assumed uniform background magnetic field 
of the surrounding medium with density Qsxt, O is the magnetic 
flux, and Qief is the background rotation rate (see also Section 
|6l). Lastly, K is the Complete Elliptic Integral of the First Kind, 
and the symbols r< and r> denote the smaller and larger of r and 
r', respectively. 

We solve equations Eqs.(fTali-(fTeb using the method of 
lines (e.g., Schiesser 1991) together with a finite volume ap- 
proach. Hereby, the spatial part of the equations are ini- 
tially discretized, and the resultant set of time-dependent or- 
dinary differential equations (OD Es) are solved with LSODE 
(iRadhakrishnan & Hindmarshlll993^) . This implici t ODE solve r 
uses the adaptive backward-difference method of iGeaJ (Il97lh . 
The smallest cell size is initially 10"^ AU and becomes 0.02 Rq 
at the highest refinement. The advection step is done using a 
second-order van-Leer TVD advection scheme dvan Leeilll977h . 
and we calculate all derivatives to second-order acc uracy on our 
nonuniform gr id (see Ciolek & Mouschoviasll 19931) . We employ 
the method of iNorman et al .I(ll980h to advect angular momen- 
tum, which avoids angular momentum diffusion and associated 
spurious ring formation. 

Our computational grid has 256 radial cells in logarithmic 
spacing. Test runs with 1024 radial cells did not show any sig- 
nificant quantitative deviations. The grid is refined to a higher 
resolution (each step by a factor of 3) whenever the column den- 
sity increases by a factor of 50. This allows us to satisfy the 
Truelove criterion dTruelove et al J 19971) . in that the Jeans length 



(I Jeansll 1 902l 1 1 9281) is resolved at all times by a minimum of 10 
cells. 

Our boundary conditions are as follows. Besides the axial 
symmetry, we have reflection symmetry at the midplane and at 
the origin. Finally, at the outer radius, we have constant-volume 
boundary conditions. For purposes of calculating derivatives we 
assume the column density to go to zero at the outer radius, 
and the magnetic field to go to its constant external value Z?,ef. 
The external medium has a low density pext and a pressure f ext- 
The assumed high temperature and low density in the external 
medium allow us to use the force-free and current-free approx- 
imation for the magnetic field above and below the cloud, as 
described earlier 



3. Magnetic braking 

We calculate ma gnetic braking using t he same technique 
presented in [ Basu & Mous choviasI (1 19941) and also used by 
iKrasnopolskv & Konigl (i2002l) . namelv an analytical approxi- 
mation to steady-state Alfven wave transport from the disk to 
an external medium. Owing to numerical complexity, a calibra- 
tion of this method with results of three-dimensional MHD wave 
propagation through a stratified compressible medium has not 
been done to date. 

In order to derive the equations for magnetic braking, we 
consider the transport of angular momentum by torsional Alfven 
waves along flux tubes from the disk to the surrounding tenuous 
medium that has the density £)ext and rotates at the angular fre- 
quency Qref . The Alfven speed in the external medium is 



Br, 



VA,ext — 



(8) 



If the timescale for the Alfven waves to reach the background 
medium far away from the disk is much smaller than the dy- 
namical evolution timescale of the disk, then the rate of angular 
momentum flux per unit radian through an annulus at radius r^^f 
of thickness drref is 



'dt 



-eext'"ref ~ ^lef) ^A,ext?"i-ef dn-gf , 



(9) 



where / denotes angular momentum per radian. The minus sign 
implies that we assume (Q - Qi-ef ) > 0, and so angular momen- 
tum is lost from the disk. In this expression we do not take into 
account any azimuthal drift of the field lines with respect to the 
neutrals, and leave that for a future study. The radius r^^f, far 
above the disk where the magnetic field has reached its uniform 
background state, is threaded by the same field line as radius r in 
the disk. Eq.(|9]l considers the angular momentum that a volume 
of gas of mass far away from the disk can take on. Angular mo- 
mentum flux is calculated by multiplying the angular momentum 
density gUn with the transport velocity VA.ext (which is constant 
far from the disk). 

We perform transformations to express Eq.(|9]i in terms of 
quantities at the disk's surface instead of the external medium. 
First, we replace the external density by the Alfven speed in 
Eq.®. Another assumption is that magnetic flux O, defined by 
Eq.dTcli. is conserved above the disk (flux freezing). Then, we 
can equate the flux in equatorial plane of the disk <J> (r) with its 
value far above the disk, where the magnetic field Z?,ef is uni- 
form. The foot point r in the disk maps to a radius rjef above, 
since the two are connected by a flux tube. We have rref > r 
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because the field lines are extending (diverging) above the disk. 
This means 



1 

d<E) - Biefiefdrref = Bc,eq'"dr. 



(10a) 
(10b) 



Using these relations, we can rewrite Eq.(|9]l to yield the angular 
momentum flux 



AJ 
d7 



27rVA,ext 



{Q. - ilref) Bc,eq'"dr. 



(11) 



Finally, taking into account the flux in both directions, above and 
below the disk (yielding a factor of 2), and switching to angular 
momentum per unit area L = SnQr- = J /rdr, we arrive at an 
expression for Nd, the torque acting on the cloud 



dL O 

-7-=N,i^ 

at ^VA,ext 



(Q-aef)Bz.e 



(12) 



At the same time, the stress-energy tensor yields the change in 
angular momentum to be equal to rB.,eqB,^,z/27r, which allows us 
to calculate the i^-component of the magnetic field at the upper 
surface of the disk as 



B. 



if.Z 



20 (Q-Qref) 

^A.ext r 



-2^^(Q-a.). 



(13) 



Together, Eqs.(fT2ll and ( fTST l yield the second term on the right 
hand side of the equation of angular momentum, Eq.dTcli. In the 
absence of magnetic braking, angular momentum is conserved 
and this term vanishes. 



4. Non-ideal MHD treatment 

We outline the derivation of the induction equation with all non- 
ideal MHD effects: ambipolar diffusion, Ohmic dissipation, and 
the Hall effect. For brevity and clarity, we only consider a sin- 
gle grain size and omit inelastic collisions. However, we include 
their effect in our calculations, and refer t he interested reader to 
the det ailed exposition in Appendix B. 1 in lKunz & MouschoviasI 
(I2OO9I) where all terms are included. 

We start with Faraday's law in cgs units, transformed to the 
frame of the neutrals: 



f..Vx(^xB-E.), 



(14) 



where B is the magnetic field, £„ is the electric field in the ref- 
erence frame of the neutrals, and c is the speed of light. In ideal 
MHD, En 5 by definition, as the conductivity is infinite and so 
afl local electric fields (in the neutral frame) are shorted by cur- 
rents immediately. This is not the case in non-ideal MHD, where 
En + 0. We therefore seek an expression for En in the general 
case. 

We take the force equations for all charged species (denoted 
by subscript s, where in our discussion s - {i, e, g'^, g"), but 
others would be possible), assuming force balance between the 
Lorentz force and collisions with neutrals. Inertial forces and 
collisions with other charged particles are neglected, as we are 
working under the assumption of a weakly-ionized plasma (with 
an ionization fraction a: 10"*^). Then we have 

= «,^,(E + ^xB)-^(v,-Vn). (15) 



The time scale for collisions between species s and neutrals is 
given by T^n- The mass density of the charged species is while 
^.5 is their charge number Note that the latter carries a minus sign 
if the charge is negative (e.g., for electrons). 

Introducing the drift velocity Wj - (Vj - Vn), we can trans- 
form Eq.lfTSll to the reference frame of the neutrals: 



= WiTjn ^^En + Wj X bj - Wj, 



(16) 



where b s B/Z? is the unit vector in the direction of the magnetic 
field, and B = |B| is the magnetic field strength. Also, we have 
introduced the cyclotron frequency (in cgs units) 



(jJs = gsB/nisC 



(17) 



where nis is the mass of the charged particle, and note that Qs s 
nisns where n , is the number density of the charged species s. 

We eliminate the term oc (wj x b) by inserting Eq.(fT6]lxb 
back into Eq.dTSll. This yields 



W( = ClJ.T, 



+ (WjTin-En X b - <y.5TsnW. 



,±), (18) 



which provides the drift velocity in terms of the electric field in 
the frame of the neutrals. 

First, we look at the component parallel to b: 



Wj,|| - WjTjn— EnJI 
D 



(19) 



We write the electric current density using the overall charge 
neuti-ahty Y^s'^s^s = 0: 



(20) 



and arrive at a generalized Ohm's law for the parallel component 
ofj 



j|| = '^n,q,co,T,„^E„^n = ^ cr^En.i 



Here, 



OS = n,qjsnlms 



(21) 



(22) 



is the conductivity due to species s. Lastly, we define cr\\-Yis '^s 
as well as t/h = l/cr|| and invert Eq.dSTb to get 



En.ll - j|| - ?7||j||- 

0^11 



(23) 



Similarly, we seek a relation between the perpendicular com- 
ponents of electric current density and electric field. Again, we 
start from Eq.lfTSll. but this time look at the perpendicular com- 
ponent: 

= f , "^'^r, ^] [K.± + W.T,nEn X b] . (24) 

Inserting this into Eq.(l20li. we find 

-Zff^Kx''- (25) 



5 



Wolf B. Dapp et al.: Disks in tlie Class phase with AD and OD 



Defining 



5. Chemistry 



.,2^2 



and cth = - > , 



Ohm's law for the perpendicular component is 

i± = cr^En,i - O-nEn X b. 



(26) 



Combining both parallel and perpendicular components, we 
can write down a generalized Ohm's law. Note that the mag- 
netic field introduces an asymmetry (and thus off-diagonal com- 
ponents) and makes a tensor expression necessary: 



J = 



cr± -CTj] ' 
(Th cr^ 

(TllJ 



E„ 



(27) 



Finally, we invert the matrix to get an expression for the electric 
field in the reference frame of the neutrals. 



En 



r]± m o\ 
7711 



(28) 



where 



cr, -Hcr- 



= 2 2 



and 7711 = — . 



The purely vertical magnetic field in our thin disk is gen- 
erated by an azimuthal current only. This means that there is 
no component of the electric current density parallel to the field 
(and neither in radial direction). Hence we are only interested in 
the perpendicular component of the electric field, and for conve- 
nience quote its expression: 



En,± = 77 + 77Hj X b, 

= mil- +{^±-m)h + miy><^- 

OD AD Hall 



(29) 



Note that the quantity 77^ contains the effects of both ambipolar 
diffusion and Ohmic dissipation. In Eq.(|29]l, the contributions of 
each of the three non-ideal MHD effects — Ohmic dissipation 
('OD'), ambipolai- diffusion ('AD'), and the Hall effect ('Hall') 
— are highlighted. We do not consider the Hall effect in the 
present work. 

Summarizing the above derivation, we can write Eq. (fT4l i as 



dt 



r dr 

Id, X Id, X 



r dr 

1 a I 



(30) 



where the relation j - c/A-nV x B has been used, and 77eff s 
T]i_c^ lAn was introduced. 



In this section, we present the chemical model used to calcu- 
late the ionization fraction, fractional abundances, and resistiv- 
ities for seven species. In our models using an MRN grain- 
size distribution there are a total of 19 species, but we omit 
the corresponding full chem ical model (which can be found in 
iKunz & Mouschoviasll2009l) for sake of clarity. In the follow- 
ing, we consider neutral matter (one helium atom per five H2 
molecules), atomic and molecular ions (such as Mg^ or HCO^), 
electrons, and grains (of uniform size, positively-charged, neu- 
tral, and negatively-charged). Multiply-charged grains are ne- 
glected, because the capture rate by a charged grain of a gas- 
phase particle with the same charge q is reduced by a factor 
exp [-q^/aksT) (ISDitzedll94Ih . where a is the particle's radius. 
A grain is thus far more likely neutralized by capturing a particle 
of opposite charge than elevated to higher charge by capturing a 
particle of the same charge. For example, .Nakano et al.. (20021) 
show that the abundance of doubly-charged grains is 5 orders of 
magnitude less than that of singly-charged ones. We also ignore 
multiply-charged ions and molecules. 

5.1. Ionization rate 

We consider four sources of ionization: 

1 . UV ionization 

2. cosmic ray ionization 

3. ionization due to radiation liberated in radioactive decay 

4. thermal ionization through collisions. 

UV ionizat ion is only important where the visual extinction 
Av < 10 (ICiolek & MouschoviasllI995h or, equivalently, the 
H2 column density Nh^ < 2 x 10^^ cm"^. The part of the 
cloud relevant to this work is much denser (in fact, everything 
within 10"* AU from the center is denser) and, to good approx- 
imation, UV ionization would not need to be considered. We 
still include its cont ribution in parametrized form , by adding 
467.64 71 cm (see lFiedler & MouschoviasllI992l) to the elec- 
tron and ion number densities. This has the effect to maintain an 
ionization fraction of ^ 3 x 10"^ in the outermost envelope of the 
core, and keep it flux-frozen, but does not affect the dynamical 
evolution of higher-density regions. 

In the higher-density regions (where 10"* cm"^ < nu^ < 
10'^ cm""*), the ionization is mainly due to cosmic rays. Their 
ionization rate is calculated by 



(cR = exp [-Eh,/ 96 g cm ^ ] 



(31) 



(lUmebavashi & Nakand[T980h . where ^0 = 5 x 10"^'^ s'' is th e 
canonical unshielded cosmic-ray ionization rate (ISpitzeilll978h . 

Beyond 77h, > 10'^ cm"^, even cosmic rays are shielded and 
cannot penetrate deeply. Here, radioactivity, mainly due to "^"K, 
still provides a background level of ionization. The ionizatio n 
rate is ^40 = 2.43 x 10"^^ s ' (e.g.. lKunz & Mouschoviasl2009l) . 
Other radionuclides (such as ^^Al) are not considered, due to 
their low abundance, their short half-life time, their low ioniza- 
tion rate, or a combination thereof. 

Finally, when the temperature reaches > 1000 K, collisions 
are energetic enough to cause thermal ionization of some atoms 
with low ionization potential (predominantly potassium). As 
the temperature rises further, collisions becomes the dominant 
source of ionization. We parametrize this as an additional source 
term in the ion equilibrium equation (see Section I5.2l i with the 
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value (see lPneuman &Mitchel]||l965h 

= 2.9 X 10"'^ cm^ s"' «H,«A» (r/1000K)'/2 

dt 

xexp(-5.03x lO'^K/r), (32) 

where n^o is the fractional abundance of the neutral atomic parti- 
cles. The fact th at potassium only makes up a; 1 /14 of all metals 
jLequeuxll 19751) has been considered in the coefficient in front. 

Beyond «Ht ^ 10'^ cm"-', at a temperature of > 1500 K, we 
assume the grains to be destroyed and the stored charges are re- 
leased into the gas. The gas beco mes highly ionized, and w e 
thus revert back to flux freezing (iPneuman & Mitchelll 11965]) . 
Note that the resistivity has already decreased significantly at 
this temperature as a consequence of thermal ionization. 

5.2. Fractional abundances 

In order to calculate the resistivity (see Section HI, we calculate 
the fractional abundance of each species using a chemical equi- 
librium network. 

The molecular ion equation describes the production of 
molecular ions (such as HCO^) by radiative ionization, their de- 
struction through charge-exchange reactions with neutral atoms 
and grains, as well as recombination reactions with electrons: 

+ nm+ng-Q'm+g- + nm+ng»Q'm+g°- (33) 

Here jS is the charge-exchange coefficient and ofdi is the dissocia- 
tive recombination rate (collisions with electrons). Their value 
and that of the other rate coefficients between species s and sf, 
as^sr, depend on temperature and grain properties, and are given 
in Appendix iBl The symbol refers to the number density of 
species s. Cosmic rays will ionize molecular hydrogen, forming 

almost instantaneously; this in turn is highly reactive and will 
strip away an electron from any (non-H2) molecule it encounters. 
For instance, CO will be transformed to HCO^ and H2. This is 
why cosmic rays act on H2 in the first term of this equation. 

Atomic ions (e.g., Na^) are produced by charge exchange 
reactions with neutral atoms, as well as thermal ionization, while 
they are destroyed by radiative recombinations with electrons, 
and by the collision with grains: 

n^+Tij^o/S + nAonHjQ'AOH, = np,+nea„ 

-H nA+Wg-ffA+g- + nA+ngoCA+g"; (34) 

where a„ is the radiative recombination rate, and the second 
term on the left-hand side represents thermal ionization (see 
SectionlSTTTl. 

The equation for positively-charged grains balances the de- 
position of charge from atomic and molecular ions with the 
capture of electrons and neutralization with negatively-charged 
grains. The charge exchange between positively-charged and 
neutral grains is in steady-state, and so their contribution appears 
on both sides of the equation and cancels: 

nni+WgOffnj+gO + ri/^+ngoap^+gO 

— «eng+a'eg+ + ng+Wg-Q'g+g-, (35) 

Similarly, negatively-charged grains form by capture of an elec- 
tron by a neutral grain, and are neutralized by charge exchange 
during collisions with molecular and atomic ions as well as 
positively-charged grains. Again, the charge exchange between 



negatively-charged and neutral grains is in steady-state, and so 
does not appear: 

neHoOO'ef" = nn,+ n(,-Q'n]+<,- -I- nA+Wg-CA+g- 



Finally, to close the system, we include equations for the total 
number of atoms and grains 

"a" + "A* = "A (37) 
«gO + «g+ + rig- - rig, (38) 

as well as for overall charge neutrality 

+ «A+ + ng+ — «g — rig — 0. (39) 

Equations (l33])-(l39b form the non-linear system to be solved 
for the fractional abundances; we do this iteratively for each 
time step and each grid point using Newton's method. The lin- 
earized matrix equation that results involves the Jacobian and 
the update vector to the abundances, and is solved directly us- 
ing a LU-decomposition package. We assume the mean mass 
of the molecular and atomic species to be - 29 and 
mA+ - 23.5 ;7Zp, respectively, where is the proton mass. Those 
values are the masses of HCO^ and the average of atomic mag- 
nesium and sodium, respectively, which are good representatives 
of the broader range of species present. 

6. Initial Conditions and Physical Units 

We assume that our initial state was reached by core 
contraction preferentially al ong magnetic field lines (e.g., 
iFiedler & MouschoviasI 19931) and rotational flattening, and pre- 
scribe initial profiles for column density and angular velocity 
given by 

So 2i2r 

2 (r) = , Q (r) = , % . (40) 

Here, R x 1,500 AU approximately equals the Jeans length at 
the core's initial central density (see below). The column den- 
sity profile is representative of t he early stage of collapse (e.g., 
lBasulll997l: Ibapp & Basull2009l) . and the angular velocity pro- 
file reflects that the specific angular momentum of any parcel is 
proportional to the enclosed mass mend- 

We assume an initial profile for B;,eq such that initially the 

normalized mass-to-flux ratio ji - In^/G S/Bj eq - 2 every- 
where, wh ich is the approximate starting point of runaway col- 
lapse (e.g.. lBasu & M ouschovias 1995J, an d also consisten t with 
observed values in molecular cloud cores (ICrutchedll99"9l) . The 
core is at rest at the beginning of our simulations, i.e., the radial 
velocity is zero. The initial state is not far from equilibrium, as 
the pressure gradient and magnetic and centrifugal forces add up 
to a; 82% of the gravitational force. Our results do not depend 
strongly on the choice of initial state, as long as gravity remains 
dominant. 

The initial central column density is set to Zq - 2 x 
10"^ g cm"^. The core has a temperature T = 10 K, and its to- 
tal mass and radius are 28.5 Mq and 0.6 pc, respectively. The 
initial central number density and magnetic field strength are 
He = 3.3 X 10"* cm"-' and B- eq ~ 200 /iG, respectively. We choose 
the external density to be next = 10^ cm"^, and the central angu- 
lar velocity Q.^ such that the cloud's edge rotates at a speed of 
0.1 km s"' pc"'. 
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Table 1. Simulation model overview. 



Model 


cigjfim 


Notes 


1 


0.019 




2 


0.038 


fiducial grain size 


3 


0.075 




4 


0.113 




5 


0.150 




6 


0.038 


no resistivity 


7 


0.038 


no magnetic braking 


8 




OD alone (as in DBIO) 


9 


MRN distribution 





" All models above start from the same initial conditions, 
comparable to those described in DBIO. For details, see 
Section|6] 



Transient adjustments occur if the simulation is started from 
an initial non-equilibrium state that is at rest. The chemistry 
calculations are quite sensitive to fluctuations in density, which 
can cause problems. We therefore let the system evolve from 
an initial state with the above-mentioned profiles and charac- 
teristics, but initially without non-ideal MHD effects. Once the 
cloud has settled into a steady infall (at a central density of 
approximately «c ~ 4 x lO*" cm"^) the full MHD equations 
Eqs.(fTab-(fTeb are solved. The state at which we switch on the 
detailed treatment of chemistry corresponds very closely to the 
initial state employed by DBIO, so that comparisons with their 
results are possible. The rotation rate by then has increased to 
1 km s ' pc~', consistent with observations of large m olecu- 
lar cloud cores dGoodman et al.ll993t[Caseni et al.l2002h . In the 
plots in this paper, we show only the region of the core within 
0.05 pc, again consistent with DBIO. Note that the nature of the 
collapse is very dynamical and happens under flux freezing to 
a very good approximation between ~ 10"* - 10'° cm""* (see 
iKunz & Mouschoviasll20 1 Oh . 

We fix the dust-to -gas ratio at 1%, consistent with obser- 
vations (ISpitzeilll978h . and keep it constant everywhere and at 
all times. This means that a different mean grain size will re- 
sult in differing total numbers of grains available, by a factor of 

" , with a total surface area oc (agr/flg,-,„in) . 

7. Results 

7.1. Collapse phase 

We ran multiple models with different single grain sizes, and 
one with a standard "MRN" distribution (named after its au- 
thors: |Mathis]^^^£L!&Nordsiec3i973^ In this case, the num- 
ber density of spherical dust grains with radii between a and 
a H- da is 

dng,tot (a) = A^MRNfl^^'^dfl. (41) 

The distribution is truncated at a lower grain radius amin = 
0.01 1 /im and an upper grain radius Omax = 0.26 fim. The coef- 
ficient AMRriis_pro2ortional^to_flw^ mass ratio in the 
cloud. iKunz & MouschoviasI (l2009h found that there is reason- 
able convergence for 5 grain size bins, so we choose that num- 
ber in order to limit computational expense. More detail can be 
found in that paper (their Section 4.2). Our models are summa- 
rized in Table [1] 

Figure [T] shows the column number density profile versus 
radius at different times for Model 2 (with Ogi- = 0.038 fim). 



Several features are identifiable via their associated breaks in the 
profile. From the outside in, those are: 

1. Prestellar infall profile with A^ oc r This corresponds to a 
volume number density profile of « oc r"^, t ypical for col- 
lapsing cores domin ated by gravity (see, e.g.. lLarso 
iDapp & Basull2009l) . The vertical hydrostatic condition man- 
dates that n oc N^. 

2. At 10 AU, a mag n etic wall dLi & McKed ^1991 
IContopoulos et al.1 119981: iTassis & M ouschovias 20oS- 
Here, the relatively well-coupled, bunched-up magnetic 
field expelled (relative to the neutrals) from the first core 
decelerates material temporarily. At smaller radii, the infall 
continues. 

3. Expansion wave profile with A^ oc r"'^^ outside th e first core. 
The power law can be motivated analytically (see lShu|[l992l 
for the spherical case). Energy conservation requires the in- 
fall speed from large distances towards a point mass (i.e., the 
first core) with mass M to scale as v,. oc y/GM/r oc r"'^^. 
At the same time the infall onto the point mass is essentially 
a steady-state process, and thus M = 2nrY.Vr - const, close 
to the border of the first core. Together, those two relations 
imply 2 oc r"'^^. 

4. First core at 1 AU. Here, the density is sufficiently high that 
the gravitational energy released in the collapse cannot es- 
cape as radiation anymore. The temperature rises, and ther- 
mal pressure gradient stabilizes the object. The first core is 
nearly in hydrostatic equilibrium, and its radial and vertical 
extent are approximately equal. 

5. Infall profile onto the second core with A^ oc r After the 
temperature in the first core has reached ^ 1 , 000 K, hydro- 
gen dissociates. This process provides a heat sink, and the 
temperature no longer increases sufficiently for thermal pres- 
sure to balance gravity. The core starts to collapse again and 
flattens at the same time. 

6. Beginning expansion wave profile with A^ oc r^'^^ outside the 
second core, for the same reasons as outside the first core. 
Once this rarefaction wave reaches the boundary of the first 
core, the material comprising the first core will fall in to a re- 
gion of stellar dimensions, unless it forms a centrifugal disk. 

7. Second core at ^ 1 Rq. After hydrogen dissociation (and 
ionization) has concluded, the temperature rises again, 
and a truly hydrostatic object, the YSO, is formed. It 
is also partly supported by electron degeneracy pressure 
dMasunaga & Inutsukall2000l) . 

Figure |2] shows the evolution of the central diffusion coeffi- 
cient T/eff versus number density for the various models (differ- 
ent single grain sizes and the MRN distribution of grain sizes; 
see Table [Hi, as we ll as the parametrized resistivity used in 
iMachida et al.l (l2006h and DBIO. At low densities, electrons and 
ions are the dominant (but not exclusive) contributors to the con- 
ductivity. Smaller grains are fairly well coupled to the magnetic 
field because of their smaller gyro-radii and smaller collisional 
cross sections, so their contribution to the effective resistivity is 
lower than that of larger grains. However, there is a competing 
effect: smaller grains have an increased capability to absorb gas- 
phase charge carriers because of their larger total surface area. 
For grains with radius Ogr = 0.075 fiva (dot-dashed line) and 
larger, the trend reverses and the conductivity increases (the re- 
sistivity decreases) with larger grain radius, as expected from an 
increased ionization fraction. 

This effect is demonstrated in Fig. [3] which shows the diffu- 
sion coefficient versus grain size at two densities (« ^10' cm"-^ 
- solid line; « ^ 10'^ cm"^ - dashed line). The symbols indicate 
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radius / AU 

Fig. 1. Column number density profile versus radius for Model 
2. The thin lines (in ascending order) are plots at the times listed 
in Table |2] Several features are identifiable via their associated 
breaks in the profile. (1) Prestellar infall profile with oc r 
(2) Magnetic wall at ^ 10 AU, where the bunched-up field lines 
decelerate material before it continues the infall. (3) Expansion 
wave profile with oc r"'^^ outside the first core. (4) First core 
at 1 AU. (5) Infall profile onto the second core with oc r 
After the first core has reached x 1 , 000 K, it starts to collapse, 
as H2 is dissociated. (6) Expansion wave profile with N cc r 
outside the second core. (7) Second core at a; 1 Rq. 



Ill 



values for which we have complete runs, while the continuous 
lines are computed by calculating the effective resistivity with 
values for temperature, column and volume density, as well as 
magnetic field averaged over the values from the available com- 



Table 2. Output times" and corresponding central number den- 
sity in the profile plots. 



Line number* 


time / lO"* yr 


«(; / cm ^ 


1 


0.0 


4.3 X 10^ 


2 


14.338 


3.3 X 10^ 


3 


19.303 


2.5 X 10** 


4 


21.063 


1.9 X lO** 


5 


21.692 


1.4 X 10^° 


6 


21.944 


1.1 X 10" 


7 


22.045 


8.6 X 10" 


8 


22.085 


6.8 X 10'2 


9 


22.105 


4.8 X 10'^ 


10 


22.116 


3.7 X 10^4 


11 


22.139 


2.8 X 10^5 


12 


22.1508 


2.1 X 10'^ 


13 


22.15137 


1.6 X 10'^ 


14 


22.15149 


1.2 X 10"* 


15 


22.15153 


9.3 X 10"* 


16 


22.15154 


7.0 X 10''^ 


17 


22.15156 


1.6 X lO^o 



" Time is counted from when chemistry and non-ideal MHD effects 

are switched on (see Section|6). 
* Lines count from bottom upwards. 



plete runs. The magnetic field strength enters the resistivity only 
though the gyro-frequency, and the estimated error in the contin- 
uous lines is typically < 1%. The diffusion coefficient takes on 
a maximum at a grain size of Ogi- ^ 0.065 ;um for low densities 
(solid line in Fig.|3]l. As described above, for grains with smaller 
radius, the effect of better coupling dominates, while at larger 
grain sizes, the associated decrease in total surface area leaves 
more free gas-phase charges, and makes for a higher conductiv- 
ity, and conversely a decreased resistivity. 

At intermediate densities the resistivity in Fig. |2] rises 
sharply, as a consequence of the grains soaking up gas-phase 
charges and decoupling from the magnetic field. Additionally, 
the conductivity drops because at high column density cosmic 
rays are shielded according to Eq. (13111. At «c ~ lO'^ cm"^ the 
only remaining source of ionization is radioactivity, which pro- 
vides a floor ionization rate (see Section 15.1b . At this stage, 
clearly distinguished by a break in the profile, the resistivity is 
almost exclusively due to grains, i.e., their collisions with neu- 
trals. Inserting Eq. dA.lb of AppendixlAlinto Eq. (l22l i. we see that 
conductivity scales as in this phase. The resistivity is the 
inverse of the conductivity, and hence it increases with larger 
grains as a consequence of their lower gyro-frequencies and 
greater collision cross sections. This happens even though there 
are fewer grains with increasing grain size, and they have a re- 
duced total surface area. The dashed line in Fig. |3] shows this 
monotonic increase in diffusion coefficient with grain radius at 
a central density of «c ~ 10'^ cm"-'. As the temperature ap- 
proaches R! 1 , 000 K, collisions become energetic enough that 
thermal ionization occurs, described by Eq.(l32li. The conductiv- 
ity recovers, and hence the resistivity decreases again. Finally, 
during the second collapse, as temperatures of « 1500 K are 
reached, grains are destroyed and all locked-up charges are re- 
leased into the gas phase. Electrons and ions flood the gas, the 
ionization fraction rises sharply, and flux-freezing is restored 
(see Section lSTb . We switch off" the chemistry calculations at this 
point, shown by the right vertical dashed line in Fig.|2]i. We point 
out that the pa r ametr i zation for only Ohmic dissipation used in 
iMachida et all ( l2006l |2007|) and DBIO (for comparison shown 
in Fig.|2]as the dotted line) yields values for the resistivity lower 
by at least a factor of 10 everywhere. 

In Model 9 (the MRN distribution), the diffusion coefficient 
is dominated by the smallest grain size present, and behaves oth- 
erwise very similar to the simulations with a single (small) con- 
stant grain size. 

Note that Fig. [2] shows data extracted from the central cell 
of a dynamical simulation. The exact value of the diffusion co- 
efficient not only depends on the grain size and the local num- 
ber density of the neutrals, but also on the number densities of 
the other species and the local values of magnetic field strength, 
temperature, and column density. As a result, the values shown 
in this plot are not necessarily representative of any other point in 
the cloud. It is therefore not helpful to parametrize the diffusion 
coefficient according to this plot and simply use an approximate 
interpolative expression as a replacement for a realistic calcula- 
tion of the microphysics and ionization balance. The procedure 
outlined in Sections |4]and|5]is the simplest method of obtaining 
accurate values for the diffusion coefficient. 

Figure |4] shows the relative contribution of ambipolar dif- 
fusion (AD) and Ohmic dissipation (OD) to the resistivity co- 
efficient near the end of the run. The coefficient is determined 
mainly by AD everywhere outside the first core at r ^ I AU. In 
fact, AD dominates OD up to a density of ~ 5 x lO'^ cm"^. 
The contribution of OD continues to increase sharply at higher 
densities and is the dominant non-ideal MHD effect within the 
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first core. Its coefficient only decreases again in the vicinity of 
the second core at ^ 0.1 AU, where thermal ionization drives up 
the conductivity. Note that the comparatively large AD diffusion 
coefficient outside 10^ AU does not cause a large flux loss 
since the dynamical time is still much smaller than the diffusion 
time associated with AD. 



10' 



lO"" 



a 10^ 



.2 10^ 



S 10' 



Sj,, = 0.019 )j,m 
a,, = 0.031! /i,m 
a,, = 0.075 lira 
ag, = 0.113 /j,m 
ag, = 0.150 i_im 
MRN distribution 
Machida ot al. 
(2006, 2007 




Fig. 2. Central diffusion coefficient //eft versus number density 
for different grain sizes, extracted from the dynamical model. 
The vertical line at the left indicates the density at which the de- 
tailed chemistry and non-ideal MHD treatment is switched on. 
Beyond ~ 10'^ cm"^ the resistivity plummets, after having 
already declined due to thermal ionization. This is where we 
switch the chemistry calculations off again, and is denoted by the 
vertical line on the right. Due to grain destruction, flux-freezing 
is restored there. 



In Figs. [5]and|6] we present the evolution of the central ion- 
ization fraction versus number density for the different models, 
and the abundances of grains and gas-phase species relative to 
molecular hydrogen. Figure |5] shows that at low densities the 
smaller grains result in a lower ionization fraction. The reason is 
that small grains are highly abundant and thus have a large sur- 
face area to which the gas-phase charges can adhere. An increas- 
ing density causes the ionization fraction to decrease with a slope 
slightly steeper than the canonical relation cc for cosmic- 
ray-dominated ionization. At intermediate densities the decrease 
in ionization fractions level off, starting at the lowest density for 
the smallest grains and at progressively higher densities for suc- 
cessively larger grains. This phenomenon is driven by the elec- 
trons adsorbing to the grains. The ions also adsorb to the grains, 
but due to a lower thermal speed than the electrons, their colli- 
sion rate is much smaller With the grains acting as a reservoir 
for charges, and cosmic-ray ionization continuing to be active, 
the decrease in ionization fraction is temporarily halted. This oc- 
curs at later stages with increasing grain size, as Fig. |6] shows, 
for the reason of the lower overall grain abundance (scaling as 
Og^), which cannot be compensated by a larger collision rate (see 

Eg. ( IB. 7b in AppendixlBl). Beyond «c ^ lO'^ cm"^, the ionization 
fraction in Fig. |5] drops precipitously for all grain sizes, as cos- 
mic rays become shielded due to high column densities. The evo- 
lution at still higher densities is determined by radioactivity, until 
finally thermal ionization kicks in at ~ 10'*' cm"-'. Potassium 
is one of the first species to be ionized, but as the temperature in- 
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Fig. 3. Diffusion coefficient versus grain size for two values of 
number density. The diffusion coefficient behaves qualitatively 
differently at different densities. While at high densities, the dif- 
fusion coefficient increases monotonically with increasing grain 
radius, it has a maximum for low densities at a grain size of 
Qgi ~ 0.065 i-im. Below this radius, the effect of better coupling 
dominates, while at larger grain sizes the decreased surface area 
of the larger grains leaves more free gas-phase charges, with as- 
sociated higher conductivity, and conversely decreased resistiv- 
ity. 
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Fig. 4. Contributions due to ambipolar diffusion (AD) and 
Ohmic dissipation (OD) to the effective resistivity in Model 2 
(grain size Og^ = 0.038 /im) at the end of the run. OD only dom- 
inates AD beyond > 10'~ cm""*, within the first core. 



creases further, grain evaporation also releases charges into the 
gas. 

Figure |7] compares the mass-to-ffux ratio (normalized to the 
value critical for collapse) for the various models at the time ther- 
mal ionization reestablishes ffux freezing, when «c k, lO'^ cm"^, 
shortly before the formation of a central protostar For compar- 
ison, we also show the mass-to-ffux ratio for a model with ffux- 
freezing throughout, and the model with only the simple pre- 
scription for Ohmic dissipation used in DBIO. The inclusion of 
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/ cm 




Fig. 6. Central fractional abundances of species. Top left: Model 1, Ogr = 0.019 /im. Middle left: Model 2, Og^ - 0.038 /im. Bottom 
left: Model 3, Ogr = 0.075 fim. The convergence of the abundances of positively- and negatively-charged grains is pushed to higher 
densities with increasing grain radius. The departure of the electron abundance from the ion abundance also happens at higher 
densities with decreasing grain surface area. The dashed vertical line at «c ~ 10* cm"^ indicates the density at which the chemistry 
and non-ideal MHD treatment are switched on, while the one at «c ~ 3 x 10'^ cm""* shows when we assume the destruction of 
grains restores flux-freezing. Top right: Model 4, Ogi = 0.113 /im. Bottom right: Model 5, flg, = 0.150 fim. Since the grain radius 
increases only by a factor of 4/3 from Model 4, the effect on the abundances is less pronounced than it was for the smaller grains. 



ambipolar diffusion and the more realistic treatment of the mi- 
crophysics lead to a further weakening of the magnetic field, 
both in the first core (by a factor of x 10'), and in the region 
outside (by a factor of ^ 2). The region with increased mass-to- 
flux ratio also extends out slightly further, also by a factor of a; 2. 
This is caused by the higher effective value of the resistivity and 
the earlier onset of its efficacy. The mass-to-flux ratio increases 
by a factor of ~ 10"^ over the course of the evolution, to x 20, 000 
times the critical value. Field strengths in main-sequence stars 
and T-Tauri stars require flux reduction of a factor a; 10^-10** 
(even if their field is dynamo-generated) if their parent molecu- 
lar clouds are assembled from subcritical fi < I atomic gas. This 
is known as the "magnetic flux problem" of star formation. Our 
results constitute substantial progress towards understanding its 
resolution. 

Figure [8] shows the evolution of the radial profile of the ver- 
tical component of the magnetic field. At low densities, in- 
creases under near flux freezing. Ambipolar diffusion is present 
and active, but is too slow to be dynamically important. Dramatic 
flux loss occurs once grains become the dominant charge carri- 
ers at n a! 10" cm"^ and the field evolution slows. This can 



be seen in the bunching-up of the thin lines that are snapshots 
at different times (at constant increments of number density, at 
times given in Table|2|. The small fluctuations near the interface 
of the first core result from differentiating rj^ff and taking a sec- 
ond derivative of B~ in the induction equation. Near the interface 
of the first core, the steep dependence of T/eft on density in the 
range 10" cm-^^ < n < 10'^ cm ^ (see Fig. O combined with 
an extreme gradient of the density, cause a jump in the diffusion 
coefficient of many orders of magnitude. 

We stress that these fluctuations are not numerical instabil- 
ities and do not grow with time. The implicit, adaptive ODE 
solver we use to solve the MHD equations with the method of 
lines (see Section|2]i is unconditionally stable. 

In Fig.|9]we present the profile of the angular velocity at vari- 
ous times. At large radii, Q oc r"' because the core evolves under 
near angular momentum conservation with the specific angular 
momentum j = Qr^ oc wignci and E oc r ' in the prestellar infall 
(cf. Fig. [TJ. Inside the expansion wave, there is a break in the 
angular ve locity profile, as expect ed from theoretical considera- 
tions (e.g.. lSaigo & Hanawall998l) . The enclosed mass in that re- 
gion is essentially the mass of the first core (the expansion wave 
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Fig. 5. Total ionization fraction versus central density for differ- 
ent grain sizes. The vertical line indicates the density at which 
the detailed chemistry and non-ideal MHD treatment is switched 
on. 
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Fig. 7. Mass-to-flux ratio yu versus radius for difi'erent grain sizes 
at the time of the formation of the second core. For comparison, 
the thin solid line shows the case with only the simplified ver- 
sion of Ohmic dissipation, as used in DBIO. In the center, the 
difference is a factor of a; 10^, but further out it it still 2x higher, 
and also shows a greater extent by about this factor. The reason 
is twofold: firstly, ambipolar diffusion is important at lower den- 
sities and, secondly, the effective resistivity is larger everywhere 
than the parametrization used in DBIO (cf. Fig.|2]l. 



itself does not contribute much), and so is approximately con- 
stant with radius. Then j oc mend ~ const, and therefore Q oc r"^. 

The radial velocity (Fig.fTOli shows the first and second cores 
very clearly. At their edges 1 AU and ^ 10"^ AU = 2 Rq, 
respectively), accretion shocks develop and the velocity drops 
precipitously. Outside the cores, in the expansion wave, the ve- 
locity follows a power law oc r"'/^. At a; 10 AU, a slight bump 
in the infall velocity hints at the magnetic wall. The fluctuations 
within the first and second cores stem from the fact that we plot 
absolute values in a log-log plot: in nearly-stable conditions as 



Class phase with AD and OD 
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Fig. 8. Evolution of B. versus radius over time. The thin lines are 
plots at the times listed in Table |2l for the fiducial grain radius 
of flgi- = 0.038 yum. 



prevalent there, the velocity can be positive or negative, but re- 
mains small. 



T 
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Fig. 9. Evolution of the angular velocity O versus radius over 
time. The thin lines are plots at the times listed in Table |2] for 
the fiducial grain radius of Ogr = 0.038 fim. The profile changes 
from being oc r ' to oc r^^ in the expansion wave, as expected 
from theoretical considerations. 



7.2. Disk formation 

The second core is a truly hydrostatic object, and very nearly 
spherical, and so the thin-disk approximation breaks down there. 
Therefore, once the central number density «c ~ 10^*' cm"^, we 
introduce a sink cell of radius a 3 Rq. Note that this is two orders 
of magnitude smaller than is commonly used (e.g. jMellon & Lil 
2008, 2009). 

In Fig. [TT]we present evidence for the formation of a cen- 
trifugal disk. The figure shows the profiles of column density. 
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Fig. 10. Evolution of the infall velocity velocity |\v| versus radius 
over time. The thin lines are plots at the times listed in Table |2] 
for the fiducial grain radius of Og, - 0.038 /im. Outside the first 
core, the magnetic wall causes a small bump. Within the first 
and second cores material is moving about inwards and outwards 
with small speeds. This causes spikes when absolute values are 
plotted on double-logarithmic axes. 



infall velocity, and the ratio of centrifugal to gravitational ac- 
celeration shortly after the introduction of the sink cell. This is 
done for Model 2, which has resistivity and the fiducial grain 
size flgr = 0.038 ;um, and for a model without resistivity. In the 
resistive model, centrifugal balance is achieved in a small re- 
gion 10 Rq) close to the center (bottom panel), while the 
flux-freezing model is braked "catastrophically", and the sup- 
port drops to minuscule values. Centrifugal balance is a neces- 
sary and suflicient condition for the formation of a bona-fide 
centrifugally-supported disk (as opposed to a larger flattened 
structure that is still infalling). At the same time all infall is 
halted there and the radial velocity plummets (middle panel), 
while in the flux-freezing model infall continues unhindered. 
After a few more months of evolution, a Toomre instability de- 
velops, and the rotationally-supported structure breaks up into a 
ring (top panel, solid line). At this point, we stop the simulation, 
because more physics would be required to follow the further 
evolution of the disk. 

DBIO also presented a model without magnetic braking 
(their Fig. 1). In that case, the entire first core turns into a large 
centrifugally-supported structure of size several AU. 

Figure [12] shows the magnetic field line topology above and 
below the disk on two scales (10 AU and 100 AU) for both the 
flux-freezing and non-ideal MHD models. The field lines are 
calculated shortly after the formation of the second core, as- 
sumi ng force-free and cur rent-free conditions above a finite thin 
disk (iMestel & Ravll985l) . The split monopole of the flux-frozen 
model (dashed lines) is created as field lines are dragged in by 
the freely falling material within the expansion wave front at 
~ 15 AU. This is replaced by a more relaxed field line structure 
in the non-ideal case (solid lines), for which the field is almost 
straight in the central region. iGaili et al.l (|2009|) presented similar 
field configurations resulting from a simplified model for Ohmic 
dissipation during the collapse. The extreme flaring of field lines 
in the flux-freezing model is a fundamental cause of the mag- 
netic braking catastrophe. 
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Fig. 11. Evidence for disk formation. Top panel: column den- 
sity profile. The system is Toomre-unstable, and breaks up into 
a ring in the resistive model (solid line), but not in the flux- 
frozen model (dotted line). Middle panel: infall profile. Infall 
is stopped in the resistive model (solid line), in innermost re- 
gion where a centrifugal disk forms, but not in the flux-frozen 
model (dotted line) Bottom panel: ratio of centrifugal over grav- 
itational acceleration. In the resistive model, centrifugal balance 
is achieved. In contrast, catastrophic magnetic braking occurs for 
the flux-frozen model, the centrifugal support drops to negligible 
values and the first core is spun down drastically. 

We can estimate the efficiency of magnetic braking by com- 
paring its instantaneous timescale with that of the dynamical 
evolution of the core. The relevant ratio is 



'^dyn 



L/Nci 

r/Vr ' 



(42) 



where L - YQ.r is the angular momentum per unit area, A^^i = 
rB^B,/2n is the torque per unit area acting on the cloud, and 
Tdyn - rjvi- is the dynamical time. Figure [T3] shows this ratio 
for catastrophic magnetic braking for our cloud after the sec- 
ond core has formed. In the region of dynamical collapse still 
characterized by a prestellar infall profile, magnetic braking is 
not very effective (its time scale being ^ 3 times that of dy- 
namic evolution) and the cloud evolves under approximate angu- 
lar momentum conservation. In the non-ideal MHD case (solid 
line), the magnetic braking time is never shorter than the dy- 
namical time, and thus catastrophic magnetic braking is avoided. 
However, the ratio comes close to unity in the expansion wave 
where a significant amount of angular momentum loss can still 
occur Within the first core however, where the magnetic field 
has been weakened by diffusive effects, magnetic braking is to- 
tally disabled, and the ratio TMB/i"dyn rises dramatically. For flux 
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Fig. 12. Magnetic field lines. The box on the left has dimensions 10 AU on each side, while the box on the right has dimensions 
100 AU. The dashed lines represent the flux-freezing model, while the solid lines show the same field lines for the model including 
non-ideal MHD effects for a grain size flg, - 0.038 yum. In both cases, the second core has just formed and is on the left axis 
midplane. The field lines straighten out significantly on small scales in the non-ideal model compared to the flux-frozen model. 



Table 3. Scaling laws in the flux-freezing case. 
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freezing (dashed line), the ratio drops below unity within the ex- 
pansion wave, and the first core marks a further sharp drop in 
the magnetic braking time to only a tenth of the dynamical time. 
Both the decline of tmb / ^dyn within the expansion wave as well 
as its sharp drop within the first core can be considered elements 
of catastrophic magnetic braking. 

In the flux-freezing case (where the magnetic field eventually 
forms essentially a split-monopole configuration, dashed line of 
Fig. [T2l i. simple analytic scaling arguments can be made to ex- 
plain the profiles. These are summarized in Table[3] The vertical 
component of the magnetic field, B., always scales like the col- 
umn density Z due to the spatially constant mass-to-flux ratio, 
and magnetic flux conservation (flux-freezing), while B, scales 
like the gravitational field. The azimuthal component of the mag- 
netic field, Bjo ^ OQr ' from the magnetic braking calculation 
(see Section[3]l. 

In the prestellar infall region under flux freezing, the mag- 
netic braking time 

Tmb = — p p ^-r oc r. (43) 



Similarly, the dynamical time in the prestellar profile scales as 

r 

Tdyn = — oc r (44) 

since iv ^ constant. The magnetic braking time and the dynam- 
ical time therefore scale with radius in the same way; in other 
words their ratio is constant. 

In the expansion wave region, the enclosed flux is dominated 
by the split-monopole, and is approximately constant (similar ar- 
guments hold for the enclosed mass). Furthermore, 2, B~, and i', 
scale differently than in the region of prestellar infall (see Table 
[3]l. Then in the expansion wave region 

Tmb oc oc r^ (45) 

Tdyn = - r^'^- (46) 

Here, the ratio of magnetic braking time and dynamical time is 
oc r'^^, as seen in Fig.[T3]between « 1 AU and x 15 AU (for the 
dashed line). This shows that the magnetic braking catastrophe 
is a fundamental property of the expansion wave in a flux-frozen 
model. In the non-ideal case, the magnetic field does not assume 
a split-monopole configuration, and neither does the mass-to- 
flux ratio remain spatially constant, and simple scaling laws are 
not applicable. 

The centrifugal radius of a mass shell, i.e., the radius at 
which the material would achieve centrifugal balance with grav- 
ity if angular momentum were conserved, can be estimated ac- 
cording to 

ref - f/GM, (47) 

where j - Qr^ is the specific angular momentum, and M is 
the mass in the central object. While this expression is strictly 
true for infall onto a point mass, we take M to be the enclosed 
mass. The result is shown in Fig. [14] This estimation is imper- 
fect within the first core, but well-justified outside it since the 
gravitational potential there looks like that of a point mass. At 
larger radii, most of the enclosed mass is in the cloud and not 
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Fig. 13. Plot of the ratio TMB/T"dyn- For both the flux-freezing 
case (dashed line) as well for non-ideal MHD (solid line), the 
magnetic braking time is reduced in the expansion wave region 
from what is was in the prestellar infall region. In the case of 
non-ideal MHD, it turns up sharply at the interface of the first 
core where diffusion has reduced the magnetic field strength sig- 
nificantly and rendered magnetic braking ineffective. In the flux 
frozen case however, the first core marks a sharp drop in mag- 
netic braking time to only a tenth of the dynamical time. 



the central mass, but by the time the expansion wave reaches 
those mass shells and rapid infall begins, approximately half of 
the e nclosed ma ss will be in the compact central object or disk 
(e.g.. [Sh^[T977l) . The solid line is the estimate based on Eq.(l47ll 
but the dotted line takes into account an estimated loss of spe- 
cific angular momentum j by a factor of 3 - 4 in the expansion 
wave, so that is reduced by a factor of 10. This estimate is sup- 
ported by data presented later in Fig. [17] and is applied to mass 
shells that are well outside the expansion wave at the time that 
we stop our simulation. The dotted line may represent an up- 
per limit to the disk radius since some more angular momentum 
may be lost by magnetic braking when the disk is forming, if 
magnetic field dissipation has not completely shut off magnetic 
braking by that time. Future calculations can settle that issue. 
Outflows will also remove some angular momentum from the 
forming disk, reducing the size further, although a counter-effect 
that we do not quantify is disk expansion due to internal angular 
momentum transfer within the disk. We assume the latter to not 
be significant in the early Class phase. The dotted line reveals 
a centrifugal radius r^f ~ 1 AU for mass shells currently at a 
distance of 200 AU from the center, and r^f ~ 10 AU for mass 
shells at 1, 500 AU. The upper axis labels the times at which the 
expansion wave (traveling at the sound speed) is estimated to 
reach these radii, and rapid infall to begin. The subsequent infall 
to the disk should take significantly less time than the expansion 
wave arrival time, since the infall is highly supersonic. Hence, 
the expansion wave arrival times can be taken to be approximate 
estimates of the central object age just as the centrifugal radii 
estimates can be taken to be upper hmits to the disk size. 

7.3. Further results 

Figure [15] shows the ratio of the radial and vertical components 
of the magnetic field, Br/B~ s^. This is a quantitative measure 
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Fig. 14. Estimated centrifugal radius of each mass shell. Beyond 
~ 10 AU, the centrifugal radius increases linearly. The solid line 
is the current centrifugal radius, while the dotted line is an es- 
timate of rcf after the expansion wave has passed and the as- 
sociated magnetic braking has reduced the specific angular mo- 
mentum (see text). The top axis shows the approximate times at 
which the expansion wave is expected to pass a certain radius. 
The symbols show the estimated centrifugal radii of 1 AU and 
10 AU. The associated times are 6 x lO"* yr and 4 x 10"* yr. 
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Fig. 15. Ratio of Br/Br sq. In the flux frozen case (dashed line), 
B, dominates everywhere within the expansion wave, indicative 
of extreme flaring and a split-monopole field configuration. With 
field diffusion active (solid line), the field lines straighten out 
quickly inside the expansion wave (cf. Fig.fTSli 



of the local bending of the magnetic field (while Fig. \T2\ is a 
visual representation). For flux freezing (dashed line), B, dom- 
inates everywhere within the expansion wave, which is indica- 
tive of the extremely-flared split-monopole field configuration. 
Foot points of magnetic flux tubes are connected to head points 
at much larger radii. This large lever arm leads to ra pid angular 
momentum loss — the magnetic braking catastrophe (iGalli et al.l 
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I2009h . In the diffusive case (solid line), the field lines straighten 
out quickly inside the expansion wave, and only a small increase 
in pitch angle occurs inside it. 

Again, we can make scaling arguments in the flux freezing 
case. The radial component of the magnetic field, Br, always 
scales like the radial gravitational field. In the prestellar region 
of the profile (outside x 100 AU) Br and B- both scale oc r ' 
and so their ratio is constant. In the expansion wave outside the 
first core (between ^ 1 AU and x 15 AU) on the other hand, 
Br oc while B- oc E oc r^'^^. Their ratio is therefore oc r^^/^. 
As before, the magnetic field does not assume a split-monopole 
configuration in the non-ideal case, and simple scaling laws are 
therefore not applicable. 




radius / AU 

Fig. 16. Ratio of diffusion time Tditf and dynamical time Tdyn at 
the end of the run. Within the expansion wave, the ratio is smaller 
than unity, and diffusion dominates over advection. 



Figure [16] shows the ratio of diff'usion time T^is and 
dynamical time Tjvn ne ar the end of our simulation run. 
IContopoulos et al.l (Il998h showed from scaling arguments that 
upon close approach to the protostar, diffusion of magnetic flux 
will eventually dominate advection. For our model, this is al- 
ready the case in the expansion wave outside the first core. The 
reason for the sudden drop in the ratio at the interface of the first 
core is that the diffusion time scale Tdifj oc Lr jriss, where L is 
the local scale length. However, not only does the density at this 
location increase very steeply, //eft itself increases very steeply 
in the relevant density range too (see Fig.|2]i. Though mitigated 
by a simultaneous decrease in the dynamical time scale (see Fig. 
[TOb . it still effects a drop of 5 orders of magnitude in the time 
scale ratio. 

In Fig. [TtI we plot the specific angular momentum versus 
enclosed mass at the initial state (dashed line) and the final snap- 
shot (solid line), as well as their ratio (dotted line). The expan- 
sion wave brings with it a maximal reduction (a factor of < 4) 
of specific angular momentum. The radius versus enclosed mass 
is plotted as the dot-dashed line. The steep rise of the latter at 
X! 4 X 10"^ Mq is a consequence of the first core dominating the 
enclosed mass, and the material within the expansion wave only 
slowly increasing it at larger radii. 




Fig. 17. Specific angular momentum versus enclosed mass at ini- 
tial state (dashed line) and final snapshot (solid line), and their 
ratio (dotted line). The expansion wave brings with it a maxi- 
mal reduction of specific angular momentum. Further insight is 
provided by the radius versus enclosed mass, plotted as the dot- 
dashed line, with radius values denoted on the second axis on the 
right. 



8. Discussion 

Our results lead us to propose the following scenario of disk 
formation and evolution in low-mass stars: 

1 . A small disk forms, initially «: 1 AU, but eventually encom- 
passes the entire first core with ^ 1 AU, since magnetic diffu- 
sion is very efficient there and inactivates magnetic braking. 

2. Figure [14] shows that the estimated centrifugal radius of a 
mass shell at radius < 50 AU lies within ^ 1 AU. That means 
that the matter can fall to this radius without hitting a rota- 
tional barrier. 

3. At the same time, Fig.[T3]shows that magnetic braking is ac- 
tive within the expansion wave, while it is dormant in the re- 
gion of dynamical infall further out. Any material within the 
expansion wave (which moves outward at the local speed of 
sound) will lose part of its angular momentum by magnetic 
braking (typically a factor of 3 - 4, cf. Fig. [TtT i. and have its 
centrifugal barrier moved further inward. This is reflected in 
the dotted line in Fig. [14] which estimates the final centrifu- 
gal radius of a mass shell. After the expansion wave passes, 
material as far as r a; 200 AU can directly accumulate onto 
the inner ^ 1 AU. Furthermore, a centrifugal radius < 10 AU 
is expected for gas currently within ^ 1500 AU, which en- 
closes about 0.4 Mq in our model. The expansion wave will 
take a! 4 X 10"* yr to reach that point. We therefore predict no 
disk larger than a 10 AU to be detectable in Class objects 
younger t han a 4 x 10" ^ yr. Th is agrees well with the obser- 
vations of iMaurv et al.l (1201 Ol) who do not find evidence for 
disks > 50 AU in Class objects. ALMA will soon aflow 
observers to test this prediction more stringently. 

4. The material in the disk will be subjected to internal 
mech anisms of angular momentum redistrib ution, e.g., the 
MRI (lBalbus&Hawlevlll991l: lBalbusll2009l) . and/or grav- 
itational torques. It is a fair assumption that a disk in 
which self-gravity is important self-regulates Xo Q ^ 1 by 
such processes (e.g.. IVorobvov & Basu l2007h . The param- 
eter Q - c^Q./ (ttGS) approximately determines whether 
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a rotating disk is unstable to fr agmentation (lToomrelll964t 
iGoldreich & Lvnden-BelJl965h . The exact critical value de- 
pends on the geometry and the equation of state, but gener- 
ally spirals and clumps w ill form if Q $ 1, while a situation 
with g > 1 is stable (e.g.. lGamnii3l2001h . 
Redistribution of angular momentum allows material in the 
inner disk to be funneled onto the star, while material fur- 
ther out gains the ex cess angular momentum and expands 
its orbit. Ifiasul d 19981) showed that internal redistribution of 
angular momentum means that the disk radius obeys 



(48) 



where Mdisk and are the disk and star mass, respectively. 
This means that the disk will expand by a factor of 10^ by 
the time that 90% of its mass has accreted onto the central 
object. 

6. By the time the central object has accreted a significant 
amount of the available material (several 10^ yr), a tenu- 
ous disk of several hu ndred AU will be present, c onsistent 
with observations (e.g.. [Andrews & Williamsll2005l) . We ex- 
pect that the disk will expand due to internal angular mo- 
mentum redistribution, but that it will occur primarily after 
the Class phase. During the Class pha se, the disk mas s 
may remain a fair fraction of the star mass (IVorobvovl201lh . 
limiting the expansion implied by Eq.(l48Tl. 

The first step o f the above scenario is broad ly consistent 
with the models of iMachida & M atsumotol (1201 ll) that have a 
small initial rotation rate. They have studied disk formation with 
their three-dimensional models that resolved the second core 
using a simplified resistivity. However, their models with rota- 
tion comparable to observations achieve a very rapidly rotat- 
ing first core, which evolves directly into a "Keplerian" (their 
words) disk before the central stellar core has been established. 
While the term "Keplerian" is normally reserved for when the 
central point mass dominates the overall gravitational potential, 
we consider that in their context it refers to centrifugal balance. 
In fact, we did obtain a similar result in DBIO for a model with 
no magnetic braking, where the first core formed with rapid ro- 
tat ion that prevented it from col lapsing further In the model 
of iMachida & Matsumotol (1201 lb a stellar core does eventu- 
ally form, but is surrounded by a Keplerian disk from the start. 
Conversely, in our models with an observationally-motivated ro- 
tation rate, the sustained effect of magnetic braking through- 
out the collapse ensures (even in the presence of Ohmic dissi- 
pation and ambipolar diffusion) that the first core can undergo 
rapid infall onto the second core. A centrifugally-supported disk 
is only built up subsequently through accretion. We start with 
a differentially-rot ating initial s tate, co nsistent with prior grav- 
itation al collapse (iBasul 1 19971) . while IMachida & Matsumotol 
(1201 ll) start with solid body rotation, which puts more angular 
momentum in the outer regions. The two possible outcomes of 
direct second collapse followed by disk formation by accretion 
versus a large centrifugally-supported disk-like first core that co- 
exists with a second core may certainly be a function of the im- 
plementation of magnetic braking and non-ideal MHD, initial 
rotation profiles, and the varieties of dimensionality and bound- 
ary conditions of different models. Future modeling can clarify 
this interesting duality of outcomes in the earliest phase of disk 
formation. 

Previous ideal and non-ideal MHD simulations of protostel- 
lar collapse focused on the for mation of large cent rifu gal disks 
of size X 100 AU. For example. lMellon & Lj(l2009l) and Li et al.l 



(1201 ll) find that the inclusion of ambipolar diffusion, Ohmic dis- 
sipation, and even the Hall term are not sufficient to avoid catas- 
trophic magnetic braking at the distances > 7 AU from the ori- 
gin that they model, for models wit h realistic initial magnetic 
field strength. IMachida et al.l (1201 l\i claim that a disk can fi- 
nally form in later stages when the envelope is depleted and 
magnetic coupling to outer moment of inertia is therefore weak. 
Krasnopolsky et al] ( 1201 Ol) invoked "anomalous" resistivity to 
weaken the effects of magnetic braking and thereby allow for 
the formation of large centrifugal disks. 

We note that there is no a priori reason why centrifugal disks 
should be large in the early accretion phase. Using an observa- 
tionally motivated magnetic-field strength and without recourse 
to "anomalous" resistivity, we have demonstrated that centrifu- 
gal disks can form during the Class stage, albeit not as large 
as a! 100 AU. Our result is made possible by our resolution 
of the innermost regions of collapse, and a sufficiently accurate 
treatment of the relevant chemical and magnetohydrodynamical 
processes. While not yet computationally feasible, well-resolved 
three-dimensional simulations that take into consideration these 
processes would represent a significant advance beyond the thin- 
disk approximation we have employed here. 

The thin-disk approximation does not allow us to follow 
jets or outflows. While those effects undoubtedly play some 
role i n re moving angular momentum from protostellar cores 
(e.g., iPud ritz & Norman 1983), they do not start to affect the 
angular momentum balance until a centrifugal disk is present. 
Angular momentum loss due to outflows after disk formation 
only strengthens our expectation that observations will not re- 
veal any disks larger than 10 AU in the early Class phase. The 
disk would be even smaller than predicted by Eqns.(l47]i and ( |48] | 
if additional angular momentum and mass were extracted by out- 
flows. 

9. Summary and conclusions 

We present results from a new axisymmetric code using the 
thin-disk approximation that calculates the collapse of rotating 
magnetized prestellar cores. While treating vertical (z) structure 
and angular momentum transport in an approximate manner, our 
code allows us to follow the evolution all the way to stellar sizes 
and near-stellar densities. We determine the abundances of seven 
different species, and consider inelastic collisions. Our effective 
resistivity includes the effect of ambipolar diffusion and Ohmic 
dissipation for five different (single) grain sizes as well as for an 
MRN distribution of grain sizes. 

By following the complex relationships between the many 
nonlinear processes at work in the formation of stars (e.g. self- 
gravity, rotation, chemistry, thermodynamics, non-ideal MHD, 
grain effects, etc.), we have been able to simultaneously ad- 
dress both the angular momentum problem and the magnetic 
flux problem of star formation. These two processes are intri- 
cately and fundamentally linked. 

We demonstrate the formation of a small centrifugally sup- 
ported disk despite the effects of magnetic braking. The "mag- 
netic braking catastrophe" is averted on small scales by substan- 
tial magnetic flux loss from the high-density region of the first 
core. This weakens the magnetic field, preventing it from spin- 
ning down the material in that region. The central mass-to-flux 
ratio increases by a factor of ~ 10"*. Shortly after the second col- 
lapse, disk formation happens very close to the central object, 
before the protostar has accreted a lot of mass {M^ < 10"^ Mq). 
This is consistent with the observational evidence of outflows at 
a very young age and the simultaneous non-detection of disks 
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> 50 AU around Class objects (iMaurvet al.ll2010l) . ALMA 
will allow observers to improve upon these constraints, by prob- 
ing for disks of size a; 10 AU. 

We propose a disk formation scenario in which centrifugal 
disks form at the earliest stage of star formation and on small 
scales. We calculate the centrifugal radius for mass shells that 
are still infalling at the end of our simulation, and estimate the 
disk to remain < 10 AU for x 4 x 10"* yr. In this case, a disk 
would not be observable around a young Class object, even 
with ALMA. However, indirect indications of the presence of 
disks — such as outflows — can still be expected. The accretion 
disk slowly grows by continued accretion but more so by internal 
angular momentum redistribution, and may eventually reach the 
size Si 100 AU observed around Class II objects. 

Acknowledgments 

The authors thank Jan Cami for providing computational fa- 
cilities to run some of the models. W.B.D. was supported 
by an NSERC Alexander Graham Bell Canada Graduate 
Scholarship. S.B. was supported by an NSERC Discovery Grant. 
M.W.K. was supported by STFC grant ST/F002505/2 during 
the early phases of this work and is currently supported by 
the National Aeronautics and Space Administration through 
Einstein Postdoctoral Fellowship Award Number PFl-120084 
issued by the Chandra X-ray Observatory Center, which is op- 
erated by the Smithsonian Astrophysical Observatory for and on 
behalf of the National Aeronautics Space Administration under 
contract NAS8-03060. 



References 

Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363 

Andrews, S. M. & Williams, J. P. 2005, ApJ, 631, 1 134 

Balbus, S. A. 2009, arXiv:0906.0854vl 

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 

Basu, S. 1997, ApJ, 485, 240 

Basu, S. 1998, ApJ, 509, 229 

Basu, S. & Mouschovias, T. Ch. 1994, ApJ, 432, 720 

Basu, S. & Mouschovias, T. Ch. 1995, ApJ, 453, 271 

Braiding, C. R. & Wardle, M. 201 1, arXiv:l 109.1370vl 

Caselli, P., Benson, P. J., Myers, P C, & Tafalla, M. 2002, ApJ, 572, 238 

Ciolek, G. E. & Basu, S. 2006, ApJ, 652, 442 

Ciolek, G. E. & Mouschovias, T. Ch. 1993, ApJ, 418, 774 

Ciolek, G. E. & Mouschovias, T. Ch. 1995, ApJ, 454, 194 

Ciolek, G. E. & Mouschovias, T. Ch. 1996, ApJ, 468, 749 

Contopoulos, 1., Ciolek, G. E., & Konigl, A. 1998, ApJ, 504, 247 

Crutcher, R. M. 1999, ApJ, 520, 706 

Dapp, W. B. & Basu, S. 2009, MNRAS, 395, 1092 

Dapp. W. B. & Basu, S. 2010, A&A. 521, L56 

Desch, S. J. & Mouschovias, T. Ch. 2001, ApJ, 550, 314 

Draine, B. T. & Sudn, B. 1987, ApJ, 320, 803 

Fiedler, R. A. & Mouschovias. T. Ch. 1992, ApJ, 391, 199 

Fiedler, R. A. & Mouschovias, T. Ch. 1993, ApJ, 415, 680 

Galli, D., Cai, M., Lizano, S., & Shu, F. H. 2009, in Revista Mexicana de 

Astronomia y Astrofisica Conference Series, Vol. 36, 143-148 
Galli, D.. Lizano. S.. Shu. F H., & Allen, A. 2006, ApJ, 647, 374 
Gammie, C. F 2001, ApJ, 553, 174 

Gear, C. W. 1971, Numerical initial value problems in ordinary differential equa- 
tions (Prentice-Hall, Engelwood Cliffs) 

Goldreich, P & Lynden-Bell, D. 1965, MNRAS, 130, 125 

Goldsmith, P F. & Arquilla, R. 1985, in Protostars and Planets U, ed. D. C. Black 
&M. S. Matthews, 137-149 

Goodman, A. A., Benson, P J., Fuller, G. A., & Myers, P C. 1993, ApJ, 406, 
528 

Hennebelle. R & Ciardi, A. 2009, A&A, 506, L29 
Hennebelle, R & Fromang, S. 2008, A&A, 477, 9 

Jeans, J. H. 1902, Philosophical Transactions of the Royal Society of London. 

Series A, Containing Papers of a Mathematical or Physical Character, 199, 1 
Jeans, J. H. 1928, Astronomy and Cosmogony (Cambridge Univ. press, 

Cambridge) 



Krasnopolsky, R. & Konigl, A. 2002, ApJ, 580, 987 
Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2010, ApJ, 716, 1541 
Kunz, M. W. & Mouschovias, T. Ch. 2009, ApJ, 693, 1895 
Kunz, M. W. & Mouschovias, T. Ch. 2010, MNRAS, 408, 322 
Larson, R. B. 1969, MNRAS, 145, 271 
Lequeux, J. 1975, A&A, 39, 257 

Li, Z.-Y, Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180 
Li, Z.-Y. & McKee, C. F 1996, ApJ, 464, 373 

Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2006, ApJ, 647, L151 
Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2007, ApJ, 670, 1198 
Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2011, PASJ, 63, 555 
Machida, M. N. & Matsumoto, T. 201 1, MNRAS, 413, 2767 
Masunaga, H. & Inutsuka, S.-i. 2000, ApJ, 531, 350 
Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 
Mauiy, A. J., Andre, P, Hennebelle, R, et al. 2010, A&A, 512, A40 
McDaniel, E. W. & Mason, E. A. 1973, Mobility and Diffusion of Ions in Gases 

(Wiley series in plasma physics) (John Wiley & Sons) 
Mellon, R. R. & Li, Z.-Y. 2008, ApJ, 681, 1356 
Mellon, R. R. & Li, Z.-Y. 2009, ApJ, 698, 922 
Mestel, L. & Ray, T. R 1985, MNRAS, 212, 275 
Mestel, L. & Spitzer, Jr., L. 1956, MNRAS, 116, 503 
Morton, S. A., Mouschovias, T. Ch., & Ciolek, G. E. 1994, ApJ, 421, 561 
Mott, N. F. & Massey, H. S. W. 1965, The theory of atomic collisions, 3rd edn. 

(Clarendon Press, Oxford) 
Mouschovias, T. Ch. 1979, ApJ, 228, 475 

Mouschovias, T. Ch. 1996, in Solar and Astrophysical Magnetohydrodynamic 

Flows, ed. K. C. Tsinganos, 505-538 
Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199 
Norman, M. L., Wilson, J. R.. & Barton. R. T. 1980, ApJ, 239, 968 
Pneuman, G. W. & Mitchell, T. P 1965, Icarus, 4, 494 
Pudritz, R. E. & Norman, C. A. 1983, ApJ, 274, 677 

Radhakrishnan, K. & Hindmarsh, A. C. 1993, Description and Use of LSODE. 

the Livermore Solver for Ordinary Differential Equations, Tech. rep., 

Lawrence Livermore National Laboratoiy 
Saigo, K. & Hanawa, T. 1998, ApJ, 493, 342 
Shu, F H. 1977, ApJ, 214, 488 

Shu, F. H. 1992, Physics of Astrophysics, Vol. II (University Science Books) 
Spitzer, Jr., L. 1941, ApJ, 93, 369 
Spitzer, Jr., L. 1948, ApJ, 107, 6 

Spitzer, Jr., L. 1978, Physical processes in the interstellar' medium (New York 

Wiley-lnterscience) 
Tassis, K. & Mouschovias, T. Ch. 2007, ApJ, 660, 388 
Toomi-e, A. 1964, ApJ, 139, 1217 
Troland, T. H. & Crutcher, R. M. 2008, ApJ. 680, 457 
Truelove, J. K., Klein, R. 1., McKee. C. F, et al. 1997, ApJ, 489, L179 
Umebayashi, T. 1983, Progress of Theoretical Physics, 69, 480 
Umebayashi, T. & Nakano, T. 1980, PASJ, 32, 405 
Umebayashi, T. & Nakano, T. 1990, MNRAS, 243, 103 
van Leer, B. 1977, Journal of Computational Physics, 23, 276 
Vorobyov, E. I. 201 1, ApJ, 729, 146 
Vorobyov, E. I. & Basu, S. 2007, MNRAS, 381, 1009 
Wardle, M. 2007, Ap&SS, 311, 35 

Watson, W. D. 1976, Reviews of Modem Physics, 48, 513 
Williams, J. R & Cieza, L. A. 201 1, ARA&A, 49, 67 



Appendix A: Collision time scales 

We compute the collision times between the diff'erent species s 
and neutrals according to the formula (see Mouschovias 1996) 



He 



(A.l) 



The quantity ks^ue is a correction factor entering the equation due 
to the fact that the gas also contains helium. The above expres- 
sion hence calculates the collision time for a charged species s 
with all neutrals. Helium contributes only a small correct ion fac- 
tor due to its low polarizab ility compared with H2 (see ISpitzeii 
1 1 978h . iMouschoviasI (1 1 996l) gives these correction factors as 



He 



1.23, ifi = /, 
1.21, if5 = e, 
1.09, if s = g+, go,oT I 



(A.2) 
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The values for the rate c on stant {crw)sH^ a re taken from 
McDaniel & Mas onI (Il973h iMott & MassevI (Il965h . and 
iCiolek & Mousch ovias respectively: 



(crw) 



'1.69X lO^^cm^ s-', ifi = /, 
1.3 X 10-''cm3 s-i, ifs^e, (A.3) 
na^iSkBT/nmH.y^ , Us = g+, go,or g^. 

The last expression assumes that the difference between the grain 
and neutral velocities is smaller than the sound speed. 



Here, is the grain mass (assumed constant), and iPgg s 1/2 is 
the probability of charge exchange between neutral and charged 
grains (both positive and negative). The probability of neutral- 
ization in Eq. ( IB. 8b is assumed to be unity. 



Appendix B: Rate Coefficients 

For convenience, we reproduce here the rate coefficients used 
in this work. They can al so be found in Appendix A of 
iKunz & MouschoviasI (120091) . 

For radiative recombination of atomic ions and electrons, 
and for the dissociative recombination of electrons and HCO^ 
ions, we respectively adopt the values dUmebavashi & Nakanol 
119901) 

an- = 2.8 X 10"'^ (300 K/Tf^^ cm^ s"' 



(B.l) 

Q-dr = 2.0 X 10"^ (300 K/r)° " cm^^ s"'. (B.2) 

For charge-exchange reacti ons between a tomic and molecular 
ions, we use the value from lWatsonl(ll976h 

/3 = 2.5x lO^'^cm^ s"'. (B.3) 

The rate coefficients pertaining to ions (both molecular and 
atomic, indicated with subscript 'i') and electrons (subscript 'e') 
on the one hand an d grains on the other are are taken from 
ISpitzeil(ll94llll948h . with refinements made bv lDraine & SutinI 
(119871) to account for the polarization of grains: 

1/2 r / _ 2 

1 + 1 



8knT 



laksT 



teg* 



1/2 



1/2 



1 + 



1 + 



2e^ 



2e^ + akuT 



8kuT 



1/2 



1 + 



1 



2e' 



2e^ + akuT 



1/2 



2aksT 
ak^T 

ak^T 



1/2 



Ve 



(B.4) 



(B.5) 



(B.6) 



(B.7) 



For the sticking prob abiUties of ele c trons or ions onto grains, we 
take the values from lUmebavashil (Il983h . = 0.6 and Pi = 
1 .0. In these equations, a is the adopted grain radius, while other 
quantities have their usual meanings. Lastly, the rate coefficients 
for charge transfer during collisions between grains are of the 
same form as the ones above, but with modified masses. 
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